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Abstract 



We describe two efficient methods of estimating the fluid permeability of 
standard models of porous media by using the statistics of continuous Brow- 
nian motion paths that initiate outside a sample and terminate on contacting 
the porous sample. The first method associates the "penetration depth" with 
a specific property of the Brownian paths, then uses the standard relation 
between penetration depth and permeability to calculate the latter. The sec- 
ond method uses Brownian paths to calculate an effective capacitance for 
the sample, then relates the capacitance, via angle-averaging theorems to the 
translational hydrodynamic friction of the sample. Finally, a result of Felder- 
hof is used to relate the latter quantity to the permeability of the sample. 
We find that the penetration depth method is highly accurate in predicting 
permeability of porous material. 
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I. INTRODUCTION 



Much theoretical effort has been expended in attempts to either estimate or bound the fluid 
permeability of a porous medium, given its average statistical properties. Techniques used 
for this include the volume-to-surface-area ratio Q, the A parameter 0, percolation ideas 
0], and concepts from the theory of Brownian motion such as mean survival distance or 
mean survival time (inverse reaction rate) []§-§|. 

One basic class of theoretical models allows one to study either packed beds or consolidated 
porous media such as sandstone. These models consist of ensembles of equal-sized imperme- 
able spherical inclusions immersed in a completely permeable medium. When used to model 
a packed bed, these spheres are non-overlapping; when used to model a consolidated porous 
medium, they are randomly located and freely overlapping. A second basic class of models 
for packed beds and/or porous media consist of periodic arrays of impermeable spherical 
inclusions immersed in a freely permeable medium. Depending on the ratio between the 
diameter of the spherical inclusions and the separation between neighboring inclusions, and 
the periodic lattice chosen, these models may consist either of overlapping or non-overlapping 
spheres. 

The use of Monte Carlo methods to solve elliptic or parabolic partial differential equations 
(PDEs) is grounded in the mathematical subject called probabilistic potential theory 
It is well known that the probability density function of a diffusing particle obeys the dif- 
fusion equation. The time-integrated, or steady-state probability density function obeys 
the Laplace equation. These facts have been greatly generalized. One can solve a large 
class of elliptic PDE's by simulating the corresponding class of (in general) biased diffusion 
processes [[<],|n]]. The diffusion Monte Carlo method is defined by the mapping of bound- 



ary value problems involving elliptic PDE's into the corresponding diffusion problems, and 
solving them by using Monte Carlo methods. This method is powerful because it can be 
implemented computationally in a massively parallel manner. 

Previous work involving one of us (J.G.) develops a class of diffusion Monte Carlo 
algorithms, the Green's function based first-passage (GFFP) methods that are very efficient 
for studying the properties of random or highly irregular two-phase media. This class of 
problems involves: 

• Interfacial surfaces that may pass very close to each other, thus creating channels that are 
narrow, but nevertheless important for understanding the bulk properties of the material. 

• Spatial intermittency: large regions with no interfacial surface; smaller regions with a 
high density of interfacial surface. Intermittency may be due e.g. to polydispersity in pore 
diameter or void size, or to a dispersion, in solution, of very large irregular inclusions, such 
as macromolecules. 

• Interfacial surfaces that may include sharp corners and edges. The difficulties these features 
cause in solving the Laplace equation are clear (In electrostatic terminology, they cause 
singularities in charge density). 



In particular, the class of GFFP algorithms, [y_l|]: 
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• Uses a purely continuum description of both sample and diffusing particle trajectory, with 
no discretization either in space or in time and hence no discretization error. This is essential 
for calculations of permeability in low-porosity materials. 

• Uses a library of Laplace Green's functions that allow a diffusing particle both to make 
large 'jumps' across regions containing no interfacial surface, and to be absorbed from a 
distance, at an absorbing boundary, in a single jump. 

• Uses 'charge-based' rather than 'voltage-based' methods of calculation, i.e., retain only the 
location of absorption events, or equivalently surface charges. This allows one to perform 
trajectory simulation using importance sampling (because each trajectory yields one surface 
charge); thus ensuring that even singularities associated with corners and edges will be 
efficiently sampled. 

This method is an adaptation and refinement of the "walk on spheres" (WOS) method of 
diffusion Monte Carlo HlOl. This method has been used since at least 1956 to treat diffusion 



problems [12fl. It was first applied in random-media problems by Zheng and Chiew |13|; 



see also 14 . Other recent studies have also generalized WOS to diffusion in two-phase 



media modeled as arrays of small homogeneous squares flTg] and cubes [|I(J;these studies 
use the simplest nontrivial Green's functions to determine first-passage probability. The 
Einstein relation may be used to obtain effective electrical conductivities ||r5|-[T7| from the 
bulk diffusion coefficient of these models. 

The present method is more efficient than WOS for diffusing particles very close to absorbing 
boundaries; as we show here, this makes it ideal for the study of low-porosity materials. It 
has been applied to solve the scalar Laplace equation, which governs both electrostatic and 
diffusion problems. 

Electrostatic applications involve the calculation of electrostatic capacitance for very gen- 
eral bodies [|i~8H 20"ll . This is an important set of problems because probabilistic potential 



theory implies fundamental relations between the capacitance and the lowest eigenvalue of 
the Laplacian in the Dirichlet problem ||19| . This in turn leads to efficient bounds or approx- 
imations for a wide variety of physical quantities in many-body systems, for example, the 
quantum mechanical scattering length in nuclear matter, or the classical rate of heat release 
in composite media. 

Diffusive applications of this class of methods involve solutions of the Smoluchowski equation, 
either in the vicinity of a single irregular object, (e.g., a ligand diffusing near a reactive site 
on a macromolecule) ||18|| , or in a system containing many absorbing objects, (e.g., systems 
with distributed trapping) ||13|| . 

Because so many physical phenomena are governed by vector Laplace equations (e.g. low- 
velocity fluid dynamics) plfl , and tensor Laplace equations (e.g. solid mechanics) (22|, it is 
natural to seek methods for extending this previous work to deal with these two classes of 
problems [|l(J. An approach to the first class is provided by the method of angle-averaging 



||23|| . This method provides an approximate relation between the translational hydrodynamic 
friction / of an arbitrary body, and the electrostatic capacitance C of a perfect conductor of 
identical size and shape. For arbitrary nonskew objects with a no-slip boundary condition, 
using angular averaging of the Oseen tensor |24]], Hubbard and Douglas (23| generalized the 



exact Stokes-Einstein result for the translational hydrodynamic friction of a sphere, 
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f sphere QtTT/Rj 



(1) 



where rj is the solvent viscosity and R is the radius of the sphere, to the translational 
hydrodynamic friction of an arbitrary impermeable body to give: 

/ = SnrjC (2) 

where C is the corresponding capacitance. This formula is exact in the "free tumbling" 
approximation, i.e., in the case that there is no coupling between the translational and 
rotational motion, and that there is sufficient noise, or disorder, in the flow to cause the 
object to occupy all possible orientations equally. But extensive tests show the method is 
accurate well beyond this set of cases. 

A central insight of the work described here is this: the angle-averaging method may be 
applied to approximate the translational friction coefficient of an "object", even if the object 
is not a connected set, e.g., if it is a collection of inclusions that constitute the impermeable 
or matrix phase of a sample of porous media. A group of well-studied models for porous 
media fit this description. The inclusions may either be randomly placed, or set in a periodic 
lattice. The former models have properties similar to many commercially important porous 
media; the latter model a set of media formed by grain consolidation. The inclusions may be 
freely overlapping or impenetrable, i.e., unable to overlap. Finally the (spherical) inclusions 
may be uniform (monodisperse) in radius, or they may be polydisperse. 

In this paper, we present two classes of accurate, computationally efficient methods of cal- 
culating permeabilities for these models and models like them. These models combine rapid 



efficient methods of simulating Brownian motion |[TT|JT^J2^] with a pair of methods for de- 
riving the permeability from the statistics of Brownian particles diffusing near a sample of 
the porous medium. The trajectories of these particles initiate outside the sample and ter- 
minate on contacting the porous matrix. The spherical collection of small inclusion spheres 
will occasionally be referred to as the porous sample in subsequent discussions. The first 
method, the "penetration depth" method, associates the fluid-dynamic penetration depth 
with a specific property of the Brownian paths, then uses the standard relation between 
penetration depth and permeability to calculate the latter. The second method, the "unit 
capacitance" method, involves using Brownian paths to calculate an effective capacitance for 
the sample, then relates the capacitance, via angle-averaging theorems, to the translational 



hydrodynamic friction of the sample. Finally, a result of Felderhof |25| is used to relate the 
latter quantity to the permeability of the sample. 

This paper is organized as follows: Section II sketches the method used here to obtain the 
electrostatic capacitance using Brownian paths. Section III develops the penetration depth 
method of calculating permeability, in which we obtain the penetration depth from Brownian 
paths, and obtain the permeability from the penetration depth. Section IV describes the unit 
capacitance method of calculating permeabilities, and applies it to the random models of 
packed beds and consolidated porous media. Section V applies both methods to determining 
the permeabilities of the random models of packed beds. Section VI applies both methods 
to determining the permeabilities of model packed beds composed of periodic arrays of 
impermeable spheres. Section VII applies both methods to determining the permeabilities 
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of polydisperse randomly packed bed. Section VIII contains our conclusions and suggestions 
for further study. 



II. CALCULATION OF THE CAPACITANCE OF A SAMPLE OF GENERAL 

SHAPE USING BROWNIAN PATHS 

In this section, we sketch the algorithm for calculating the electrostatic capacitance of an 
irregularly shaped conducting body by using Brownian trajectories that initiate on a spher- 
ical launch surface surrounding the conducting body, and terminate on contact with that 
body. We focus this derivation in several ways in order to apply it to samples of porous 
media. If a packed bed model is either randomly packed or densely packed, any smooth 
sample boundary will intersect some of the spherical inclusions that constitute the sample; 
we provide two methods for dealing with this problem in section IV, the "effective radius 
method" and the "sharp boundary method". Difficulties will occur if the sample is chosen 
either too large or too small; we give an analysis that yields the acceptable qualitive range 
of sample sizes in section IV. For periodic models of porous media, a non-round, e.g. cubic, 
sample (and thus a non-round launch surface) might be more appropriate. The derivation 
provided here is readily generalized to allow a cubic launch surface; in forthcoming work we 
explore this possibility numerically. 

The electrostatic potential 0(r), at position r, in the presence of a conducting object, and 
the probability density p(r) in the associated Smoluchowski diffusion problem, i.e., the 
probability associated with finding a diffusing particle at position r, in the presence of an 
absorbing object of size and shape identical to the conductor, are related by: <p(r) = 1— p(r). 
This equation is derived by noting that its LHS and RHS obey the Laplace equation with 
identical boundary conditions. 

Thus the capacitance C of a conducting object is given by: 

C = -(An)- 1 I da ■ V0(r) = (An)- 1 [ da ■ Vp(r). (3) 
Jn Jn 

where fl is the surface of the object. 

The diffusion-controlled reaction rate k is defined as |I8| , |26| , |2T 

k = D jda- Vp(r) (4) 



where the integral is evaluated on any closed surface containing the inner boundary Vt 
(because of Green's theorem any such surface may be used.). 

The relation between the diffusion-controlled reaction rate, in the Smoluchowski problem, of 
particles absorbed upon contact with an object, and the electrical capacitance of a conductor 
with identical size and shape to that object is thus given by: 

k = AnDC (5) 
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We define Poo(r;£l) to be the probability that the diffusing particle started at r will be 
absorbed on Q and 700 (r;0) the probability that the diffusing particle started at r will go 
to infinity, i.e., that it will not be absorbed in finite time. Thus: 

7oo(r;fi) = l-/3oo(r;fi). (6) 



7oo(r; n) can be shown to satisfy the Laplace equation |^7|]. If the obvious boundary condi 



tions, 700(00; fi) = 1 and 7oo(^*; fi) = for r on Vt are considered, it is clear that 

P(r) = 7oo(r;fi). (7) 

Using Eqs. ^| and |7|, Eq. [| becomes: 

k = —D <j> do ■ V/3 00 (r; fi) (8) 

= -D j>r 2 sm.9d9d(j)— /Mrjfi), (9) 

= -i*D—ffsiied0d<i>/3 oo (r-,n), (10) 



Defining 



we obtain from Eq. |l0| 



P(r]Q) = — ismededtPnfcn), (11) 

All J 



d 

— 47ir 2 D—j3(r; Q) = k. (12) 



For any spherical surface r = b which contains the boundary Q, integrating Eq. |l2| from b 
to r with respect to r and using the boundary condition /3(oo; Q) = 0, gives an expression 
for the reaction rate k: 

K = 4irDb[3(b;n) (13) 

Comparing Eqs. || and [13] yields an important formula for the electrostatic capacitance 

C = bp{b; Q) (14) 

= bp (15) 
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where (3 is the fraction of diffusing particles started at a random, i.e., angle-averaged, position 
on the launch sphere that are absorbed on the target. The unit capacitance method, as 
described here, uses this relation to determine the capacitance of a packed bed by performing 
Monte Carlo simulations for the quantity (3. Note that no assumption is made in the above 
that the "object" being studied is a connected set; it may be taken to be the set of inclusions 
that constitute the matrix phase of a sample of porous media. 

We calculate the quantity (3 by performing simulations as follows: each random diffusion path 
is constructed as a series of first-passage propagation jumps, each from the present position 
of the particle to a new position on a first-passage surface drawn around the present position. 
The new position is sampled from a first-passage position distribution function f(x) which 
gives the probability density associated with finding a diffusing particle leaving the present 
position and first contacting the first-passage surface at point x. 

This simplest first-passage surface is the first-passage sphere, i.e., a large sphere centered 
around the present position of the diffusing particle, that does not intersect any of the 
inclusions. The first-passage position density on this surface is uniform, i.e., its distribution 
function is trivial. Using only this first-passage surface yields a trivial case of the GFFP 
method, namely the WOS method. In this method, whenever a diffusing particle gets close 
enough (within a fixed distance 5h) to one of the inclusions, it is taken to be absorbed. 

We show here that for the problems under study it is more efficient to use more complex 
first-passage surfaces. We let the first-passage sphere intersect the nearest inclusion, and 
grow as large as possible, provided that it: 

• not intersect the next-nearest inclusion 

• not intersect any corners or edges of the nearest inclusion 

The resulting first-passage surface is the portion of the first-passage sphere not contained in 
any inclusion. Its surface consists of part of the first-passage sphere, and part of the surface of 
the inclusion. The probability distribution function f(x,P), for a diffusing particle starting 
at the center of the first-passage sphere and making first-passage at the point x is in general 
quite nontrivial. (Here P is the set of geometric parameters that characterize the first- 
passage surface. Also, we assume for simplicity, that symmetry allows specification of the 
first-passage position a; by a single position parameter, the polar angle 9.) Thus, (P, 9) is 
the full set of parameters that the first-passage distribution function depends on. 

For a wide class of inclusions, the resulting first-passage surfaces are 'locally simple' closed 
geometrical surfaces. These surfaces are those for which we can tabulate, as a function of 
several parameters, the corresponding first-passage distribution function. To obtain this 
function, we: 

• Calculate the gradient, normal to the absorbing surface, of the appropriate Laplace Green's 
function for a grid of parameter and field values. 

• Calculate the indefinite integral of this gradient as a function of cos#. Normalize the 
integral to be unity: 

J(P,cos0) = l (16) 

for cos 9=1. 
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• Invert this relation for each set of parameter values P to obtain cos 6 as a function of P 
and /. The first-passage position, i.e., cos^, can then be importance-sampled by choosing a 
random number a G [0, 1], setting I = a, and interpolating this relation. 

If the chosen first-passage position is on the portion of inclusion surface included in the first- 
passage surface, the particle is absorbed. Otherwise, it will reach a point on the rest of the 
first-passage surface. From there it makes its next jump. The diffusing particle continues 
until it is absorbed in the absorbing target or goes to infinity. 

As shown in Fig. [I], the GFFP method is more efficient than the WOS method whenever a 
small distance 5h must be used; this is true because the average number of moves required 
for absorption in the WOS method increases rapidly as 5h decreases |10|. Thus, it is more 
efficient for studying porous media either at low porosities or high polydispersity. 

The set of tabulated first-passage distribution functions is built up by a process of 'bootstrap 
diffusion Monte Carlo', in which simple distribution functions are used to accelerate the 
simulations used to calculate more complex distribution functions, and so on. In practice, 
the set of distribution functions is limited mainly by our inability to store in cache memory 
a tabulated distribution function depending on more than three parameters. In the present 
study the inclusions are asymmetric lenses (in the sharp-boundary method) and spheres. 
These are readily treated by the GFFP method HlT| . |28| . 



III. THE PENETRATION DEPTH METHOD OF CALCULATING 
PERMEABILITY FROM THE ABSORPTION POINT STATISTICS OF 

BROWNIAN TRAJECTORIES 



In this section, we describe the 'penetration depth' (PD) method for calculating permeability 
of an absorbing sample. First we relate the fluid-dynamic concept of "penetration depth" to 
a property of the paths of Brownian particles that diffuse from outside the sample of porous 
media and terminate on contact with it. Then we use the standard relation between the 
permeability of a sample and its penetration depth to calculate the latter. 

A simple model calculation allows one to relate the permeability of a porous medium to the 



effective penetration depth of fluid into the medium fl24|| . Consider a homogeneous porous 
medium in a half-space with constant permeability k. Fluid permeates this medium and has 
a constant velocity V within it. Hz is the distance between the macroscopic fluid-medium 



interface and a point of the porous medium, the Debijf-Brinkman equation |29| 

1 
k 



VP = r/V 2 V - ^V. (17) 



where V and P are macroscopic velocity and pressure respectively, 
reduces to the following, 



' = 1 1 ' (18) 
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with the solution 



V(z) = V exp(-z/Vk) 



(19) 



This result shows that \/k measures the distance that the flow effectively penetrates into 
the porous medium. We will call \/k the 'penetration depth'. 

We can use this result to measure the permeability of a sample by identifying the penetra- 
tion depth with the difference Z, between the average radial position at which the diffusing 
particles are absorbed and the actual sample radius, thus yielding the approximate relation: 

k = I 2 (20) 



The penetration depth defined here for diffusing particles is different from the mean linear 
survival distance: the latter is the average distance from the random starting point in the 
void phase to the absorption point; the former is the average of the radial component of this 
distance. It is possible that the mean linear survival distance can give an upper bound on 



permeability analogous to the Torquato-Kim upper bound |14[ based on the mean survival 
time. 



IV. THE UNIT CAPACITANCE METHOD OF CALCULATING PERMEABILITY 

OF RANDOM SYSTEMS 

In this section, we describe an algorithm, the unit capacitance (UC) method, for calculating 
the permeability of a sample of random media. We calculate the capacitance of a sample 
(as described in Section II), use Eq. |[ from angle-averaging to give the translational hy- 
drodynamic friction in terms of the capacitance, and finally use a relation, first published 
by Felderhof, between the translational hydrodynamic friction and the permeability to de- 
termine the latter. 

Assuming that a given spherical sample is much larger than either the average distance 
between spherical inclusions or the correlation length associated with the statistics of the 
packed bed, it can be modeled as a homogeneous porous sphere with the appropriate porosity. 
One can solve the linear Stokes equation P^ , [2^| for the translational frictional coefficient, 
/, of such a sphere by making the assumption that the excess pressure is linear in the fluid 
velocity, and by using the symmetry of the problem. The result is: 

/ = 6n V RG (a){l + A Go(a)} -i (21) 

where rj is the fluid viscosity and the function G (a) is given by: 

G (a) = 1 - -tanho". (22) 
a 
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Here a is the dimensionless quantity defined by 

R 

a = 7P (23) 

where R is the porous sphere radius and k is the permeability. 

Eliminating the translational frictional coefficient between Eqs. |] and [21], one finds a relation 
between the capacitance and the permeability 

^ = G (a){l + JjGoC*)}" 1 . (24) 
Obtaining C/R (unit capacitance) from simulation allows us to use Eq. |1 to calculate a, 



and thus obtain the desired permeability estimate from Eq. |23. 

An important finding of the present work is that, even though neither the graph of C/R 
vs. R, nor that of a vs. C/R contains a flat region, the composition of the two, giving our 
permeability estimate as a function of R, does contain a flat "plateau" region. For example, 
Fig. H shows that for a broad range of porosities there is a plateau range of sample radii over 
which the permeability estimate does not vary. The lower and upper limits of this region 
vary systematically. 

In monodisperse homogeneous porous media models, the permeability is dependent only on 
the porosity (defined as the ratio of void space volume to total volume) and the inclusion 
sphere radius. By dimensional analysis, it can be expressed in the following functional form: 

k = /(0i)a 2 (25) 



where (fix is porosity and a is inclusion sphere radius. The permeability estimates produced 
by the UC method will, in addition, depend on the parameter R/a. However, once this 
parameter is set larger than both the correlation length of the medium (very small in random 
media) and the average distance between inclusions, this dependency should become quite 
weak. 



In Eq. [24], for a given porosity and inclusion sphere radius (i.e. the permeability is constant), 
if we increase the sample radius, the unit capacitance goes to unity as cr goes to infinity. 
For each porosity, there is a range of sample radii which are much bigger than the average 
distance between spherical inclusions, but for which the ratio C/R will be far enough from 
unity to permit us to interpolate the corresponding a value. We choose a radius that is in 
this range. 

When C/R is close to unity, that is, for very low porosities, the a interpolation becomes 
unstable: a small change in unit capacitance causes a large change in a. This makes the 
permeability estimate in the unit capacitance method difficult; thus the penetration depth 
method will prove to be more reliable at very low porosities. 

A fundamental problem of this study is that of defining the sample boundary. To see this, 
note the following paradox: we study the bulk property (permeability) of a medium using 
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diffusing particles that start outside a finite sample and are absorbed on its surface. These 
particles will in general be absorbed in a surface boundary layer that grows thinner as the 
porosity decreases. Thus, we face the problem of freeing our simulation from surface artifacts 
so that it can determine this bulk property, when the method is based on events that occur 
near the surface! In this paper, we show that this is indeed quite possible. 

We generate samples of the random media models that we study using two different methods: 
the "effective-radius" method, and the "sharp-boundary" method. We describe each method 
in turn. 

In the "effective- radius" method, if randomly overlapping, i.e., when Poisson statistics are 
used, we place inclusion sphere centers at random in the sample sphere of radius R. We 
continue to place centers until the density of sphere centers approaches p, where p and the 
porosity <p\ being modeled are related by: 

01 = exp[-(47r/3)i? 3 p] (26) 

If random non-overlapping statistics are used, we place centers sequentially at randomly 
chosen positions, requiring only that each center placed be far enough from the centers 
already placed that the corresponding inclusions not overlap. The resulting distribution is 
not identical to true random non-overlapping statistics, i.e., to what liquid-state theorists call 
the "hard-sphere" distribution, but studies of these two distributions show that differences 
in bulk properties emerge only at very low porosities. Even though an inclusion extends in 
part outside the sample sphere of radius R, it still contributes both to the density of centers 
p and to the simulation. 

Because of the lack of the contribution of the inclusion spheres in the region r e [R, R + 2a] , 
the region r G [R — a,R + a] will have local porosity less than the desired bulk porosity. This 
is an important concern: especially at lower porosities, the diffusing particles will seldom 
sample any part of the sample except for the boundary layer. Because of this effect, an 
effective sample radius, which is expected to be less than R, is used. For the non-overlapping 
case, we use for an effective radius the radius of a sphere, which, if it contained the same 
number of surviving inclusion centers as our sample, would have as its average porosity the 
same as the bulk porosity of our sample, i.e., the local porosity far away from the sample 
boundary. (In previous research p0|| , we used this choice of effective radius for both the 
overlapping and non-overlapping case.) For the overlapping case, the same procedure is 
used; the set of effective sample radii is then rescaled so that the effective radius at the 
critical porosity (0^=0.03) makes the calculated unit capacitance equal to 1. (Here by the 
term 'critical porosity', we mean the porosity below which fluid ceases to flow through the 
sample, i.e., the percolation threshold.) 

In the "sharp-boundary" method, instead of using an ad hoc effective radius, the porous 
media sample is constructed as follows: we place the centers of inclusion spheres into a 
large sphere of radius (R + a) according to the chosen statistics, for a given porosity. We 
then define the actual sample by drawing a sphere of radius R and allowing it to freely 
intersect inclusions already placed. The sample is then defined to be all of the void phase, all 
inclusions, and all fragments of inclusions, that are contained in this sample sphere of radius 
R. With this sharp boundary, the porosity in the actual sample boundary is maintained 
uniform up to the boundary. 



11 



The sampling of periodic grain consolidation models is also done according to the sharp 
boundary method: we take a spherical sample of radius R after making periodic grain 
models in a cube which is slightly bigger than 2R x 2R x 2R. For these models, the result 
will depend upon the position of the sample center; thus we average the results over ten 
samples, each with a randomly chosen center. 



V. STUDIES OF RANDOMLY PACKED BEDS 



In this section, we apply the two permeability estimation methods to the calculation of 
permeability of model packed beds, composed of non-overlapping, randomly placed, impen- 
etrable spherical inclusions, and model porous media, composed of randomly placed, freely 
overlapping, impenetrable spherical inclusions. We compare our results with the available 
numerical solutions of the Stokes equation, and also with a number of theoretical bounds 
and estimates from the literature. 

Our simulations show that, for the overlapping sphere model studied here, the two perme- 
ability estimation methods, the penetration depth method and the unit capacitance method, 
agree well within random error, estimated, in this the relative standard deviation 

of the results. For the nonoverlapping sphere model, the two methods deviate as poros- 
ity decreases. For both models, we obtain 10 packing bed samples and average the ten 
permeability estimates. 

As already discussed, the graph of permeability shows a very substantial plateau as a func- 
tion of sample radius, i.e., a region of sample radius over which the permeability estimate 
shows almost no variation. Also, with the PD method the same property of sampling-size- 
independent estimation is observed (see Fig. [|). 

Our estimation results for the randomly non-overlapping model are compared with some 
other data sets in Fig. §. The other data sets used, numbered a) to d), are: 



a) The Stokes' law f21fl for a dilute bed of spheres is 



9(1 -<h) 



(27) 



where a is the inclusion sphere radius. 



b) The Happel-Brenner approximation [21 

,3 



is. 



k 



■7 + h 5 - 37 e 



3 + 2 7 £ 



(28) 



where a is the inclusion radius, and 7 3 = (1 — 0i). 



c) The upper bound of Torquato-Kim [[H 

2 



is: 



k 



9t (1 



(29) 
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where D is the diffusion constant and r is the average survival time. The data for r is taken 



from Ref. [Tj. 



d) The empirical Kozeny-Carmen relation |32|.|3"3|| in the general case is 



k = Tm (T^W" 2 < 30) 

The Stokes' law, which is the simplest approximation for the dilute bed, gives estimates that 
are greater than all other data sets. The Happel-Brenner approximation and the Torquato- 
Kim upper bound are above our data. Our simulation data with the effective sample radius is 
above the Torquato-Kim upper bound. Our simulation data with the sharp sample boundary 
lies between the Happel-Brenner approximation and the empirical Kozeny-Carmen relation. 
Because our simulation data with the effective radius method violates the Torquato-Kim 
upper bound, the sharp boundary method is clearly to be preferred for this case. 

Our estimation results for the overlapping case R = 15.0 in Fig. are compared with some 
other data. The estimation is not applied at very low porosities (0i < 0.1) and very high 
porosities (0i > 0.9). 

The other data sets used, numbered e) to g), are as follows: 

e) The Torquato-Kim upper bound on permeability |3lH for an overlapping sphere bed is, 

k < -y (31) 

where r is the average survival time and F is the formation factor. Data for r and F are 



obtained from Refs. 14 and iBlL 



f) The Doi upper bound @] on permeability for this system is 

k = \a 2 e^ f 1 dx{±- + x(l - xf) exp(- 7 (^ - ^x 3 )), (32) 

where 7 is given in terms of the porosity by <p\ = exp(~7). 

g) The lattice-Boltzmann data p5j| , which agrees well with the Kozeny-Carmen equation, 

k = 0! 3 /6s 2 . (33) 

The lattice-Boltzmann method is a variant of the lattice-gas automata method. The specific 
surface area s is given analytically by s = 4^a, 2 p4>i where p is the number of sphere centers 



per unit volume. At low porosity, finite size effects occur [|36 . 

Simulation data with the sharp boundary sampling method are below the two upper bounds, 
the Doi upper bound and the Torquato-Kim upper bound. Note that the lattice Boltzmann 
result at high porosities and the actual data for monosized sphere beds at low porosities 
agree well with our simulations. However, the PD method overestimates permeabilities at 
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very low porosities because this method does not take into account the phenomenon of 
percolation. Our data from the effective radius method is quite good for the overlapping 
case. 

At extremely low porosities, estimation instability appears; when the porosity is very close 
to unity, simulation may give a sample unit capacitance slightly above 1.0 due to the Monte 
Carlo fluctuations. The unit capacitance data which exceed unit value are discarded. The 
simulations described here are large enough that this effect is quite small except at some 
low porosities of big porous samples like R = 50.0. 



VI. STUDIES OF PERIODIC GRAIN CONSOLIDATION MODELS OF POROUS 

MEDIA 



In this section, with the sharp boundary sampling method, we apply the two methods of 
determining permeability that we have developed to the permeability calculation of model 
packed beds, composed of non-overlapping, periodic, impermeable spherical inclusions, and 
model porous media, composed of periodic overlapping, impermeable spherical inclusions. 
The difference between the two types of models consists only in the ratio of the sphere 
diameter to the separation between the spherical inclusions. We compare our results with 



the solution of the Stokes equation, using fluid dynamics codes, by Larson and Higdon |3T 
for simple cubic (SC), body-centered cubic (BCC), and face-centered cubic (FCC) lattices. 

In all cases, 10 packed bed samples are used, and the results are averaged. Figures || [7] 



and H compare results from the two methods with the Larson-Higdon |37[ calculations for 
sample radius R = 15.0. At low porosities, the standard deviation for the unit capacitance 
estimation is higher because of the estimation instability. With the sampling radius R = 
15.0, two estimation results, UC and PD, show some deviations. The PD estimate is better 
than UC estimate. However, with R = 50.0, the deviation between the two methods becomes 
smaller. At low porosities, due to the estimation instability of UC, no UC data are shown. 
It seems that for grain consolidation models UC estimation needs larger sampling due to 
the cubic lattice structure. 

Based on these results, we increase the sampling size to R = 50.0. As already discussed, we 
cannot use the UC method here. The results from the PD method, Figs. §, [l(| and [TT] are 
in excellent agreement with the Larson-Higdon results. 

The same property of sampling-size-independent estimation as the random media is observed 



in Fig. O. It should be noted that, for high porosities, the necessary sample radius can be 



as large as R = 50.0. 



VII. STUDIES OF POLYDISPERSE RANDOMLY POROUS MEDIA 



In this section, we apply the two new methods with the sharp boundary sampling method 
to the calculation of the permeability of model porous media, composed of polydisperse 
overlapping, randomly placed, impermeable spherical inclusions. The inclusion sphere radii 
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are chosen at random from the values {1.5, 3.5, 5.5, 7.5}. We compare our results with the 
available numerical solutions of the Stokes equation fl36l] . 

Figure [L^ shows that the two sets of results agree well. At very low porosities, our methods 
give permeabilities greater than those given by the Stokes equation because we have not yet 
incorporated the phenomenon of percolation. This will become important when we perform 
more extensive study of polydisperse models; the percolation threshold seem to be reached 
at lower porosities in such systems. 

VIII. CONCLUSIONS AND SUGGESTIONS FOR FURTHER STUDY 

Our methods give very good results for all models of porous media tested. The penetration 
length method is better at low porosities; the unit capacitance method shows high standard 
deviations at low porosities due to the steepness of the unit capacitance curve. It is very un- 
fortunate that so little high quality simulation data exists even for the simple models studied 
here. We note that these models have been standards for theoretical study for decades. Our 
method will predict permeabilities for a large class of homogeneous and isotropic porous 
media, in the medium and high porosity regimes (0.1 < <pi < 0.9). 

The unit capacitance estimation discrepancy with sampling radius R = 15.0 may be due 
to the importance, for periodic models, of using a launch surface with the geometry of the 
sample, e.g., a cubic launch surface for SC samples. We have laid the basis for this step 
here; it will require further study. 

An important point of this method is that it is very fast compared with other methods. 
Our computations were in parallel Fortran, using MPI, on a PC cluster (32 node, 233 Mhz 
machine). Using 10 nodes, one porous sample per node, one million diffusing particles 
for each porous sample and porosity, and 10 porosities, it takes about 10 hours per set 
of calculations. Using e.g. boundary element/finite element methods to solve the Stokes 
equation in a sample of porous media can require the same amount of time to do this set of 
calculations for a single value of porosity (We note that the latter methods are, at present, 
more efficient than these methods developed here for obtaining the detailed flow field in a 
sample; this may be important in particular applications.). 

ACKNOWLEDGMENTS 

We give special thanks to Nicos Martys, Eduardo Glandt and Dan Rothman for their useful 
information and comments. We thank Joe Hubbard and Jack Douglas for useful discussions. 
We are especially grateful to Nicos Martys for sharing with us the raw simulation data used 
in reference ||36|| . 

Also, we wish to acknowledge NCSA, especially the Condensed Matter Physics section of the 
Applications Group, for access to their computational resources and collaborative support. 
In addition, we thank NCSA's educational programs for training support to use NCSA 
computational resources. 



15 



REFERENCES 



[1] A.E. Scheidegger, The Physics of Flow Through Porous Media (University of Toronto 
Press, Toronto, 1974) 

[2] D.L. Johnson, J. Koplik, and L.M. Schwartz, "New pore-size parameter characterizing 

transport in porous media," Phys. Rev. Lett. 57, 2564 (1986) 
[3] A.J. Katz and A.H. Thompson, "Quantitative prediction of permeability in porous rock," 

Phys. Rev. B 34, 8179 (1986) 
[4] A.E. Scheidegger, "Statistical hydrodynamics on porous media," J. of Appl. Phys. 25, 

994(1954) 

[5] S. Torquato, "Relationship between permeability and diffusion-controlled trapping con- 
stant of porous media," Phys. Rev. Lett. 64, 2644 (1990) 

[6] M. Avellaneda and S. Torquato, "Rigorous link between fluid permeability, electrical 
conductivity and relaxation times for transport in porous media," Phys. Fluids A 3, 
2529 (1991). 

[7] M. Friedlin, Functional Integration and Partial Differential Equations (Princeton Univ. 
Press, N.Y., 1985) 

[8] K.L. Chung and Z.X. Zhao, From Brownian Motion to Schrddinger's Equation, 

(Springer- Verlag, Berlin, 1995) 
[9] A.-S. Sznitman, Brownian Motion, Obstacles, and Random Media (Springer- Verlag, 

NY, 1998) 

[10] K.K. Sabelfeld, Monte Carlo Methods in Boundary Value Problems (Springer- Verlag, 
Berlin, 1991) 

[11] J. A. Given, J.B. Hubbard and J.F. Douglas, "A first-passage algorithm for the hydro- 
dynamic friction and diffusion-limited reaction rate of macromolecule", J. Chem. Phys. 
106(9), 3761 (1997) 

[12] M.E. Mtiller, "Some continuous Monte Carlo methods for the Dirichlet problem," Ann 

Math Statistics, 27(3), 569 (1956) 
[13] L. H. Zheng and Y. C. Chiew, "Computer simulation of diffusion-controlled reactions in 

dispersions of spherical sink," J. Chem. Phys. 90, 322 (1989) 
[14] S. Torquato, and I. C. Kim, "Cross-property relations for momentum and diffusional 

transport in porous media," J. Appl. Phys. 72(7), 2612(1992) 
[15] R.A. Siegel and R. Langer, "A new Monte Carlo approach to diffusion in constricted 

porous geometries," J. Colloid and Interface Science 109, 426 (1986) 
[16] S. Torquato, I.C. Kim, and D. Cule, "Effective conductivity, dielectric constant, and 

diffusion coefficient of digitized composite media via first-passage-time equations," J. 

Appl. Phys. 65, 1560 (1999) 
[17] L.M. Schwartz and J.R. Banavar, "Transport properties of disordered continuum sys- 
tems," Phys. Rev. B 39, 11965 (1989) 
[18] B.A. Luty, J. A. McCammon, and H.-X. Zhou, "Diffusive reaction rates from Brownian 

dynamics simulations: replacing the outer cutoff surface by an analytical treatment," J. 

Chem. Phys. 97, 5682 (1992) 
[19] J.F. Douglas and E. Garboczi, "Intrinsic viscosity and the polarizability of particles 

having a wide range of shapes," Adv. Chem. Phys. 91, 85 (1995); E. Garboczi and J.F. 

Douglas, "Intrinsic conductivity of objects having arbitrary shape and conductivity," 

Phys. Rev. E 53, 6169 (1996) 



16 



[20] H.-X. Zhou, A. Szabo, J.F. Douglas, and J.B. Hubbard, "A Brownian dynamics algo- 
rithm for calculating the hydrodynamic friction and the electrostatic capacitance of an 
arbitrarily shaped object," J. Chem. Phys. 100(5), 3821 (1994) 

[21] J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics ( Noordhoff Interna- 
tional Publishing Leyden, Netherlands, 1973) 

[22] K.K. Sabelfeld, "Integral and probabilistic representations for systems of elliptic equa- 
tions," Math. Comput. Modelling, 23, 111 (1996) 

[23] J.B. Hubbard and J.F. Douglas, "Hydrodynamic friction of arbitrarily shaped Brownian 
particles," Phys. Rev. E 47, R2983 (1993) 

[24] F.W. Wiegel, Fluid Flow through Porous Macromolecular Systems (Springer- Verlag, 
Berlin Heidelberg New York, 1980) 

[25] B. U. Felderhof, "Frictional properties of dilute polymer solutions; III. translational- 
frictional coefficient," Physica A 80, 63 (1975). 

[26] H.-X. Zhou, "Calculation of translational friction and intrinsic viscosity. I. general for- 
mulation for arbitrarily shaped particles," Biophysical Journal, 69, 2286 (1990) 

[27] H.-X. Zhou, "On the calculation of diffusive reaction rates using Brownian dynamics 
simulations," J. Chem. Phys. 92(5), 3092 (1990) 

[28] H.M. MacDonald, "The electrical distribution on a conductor bounded by two spherical 
surfaces cutting at any angle," Proc. London Math. Soc. 26, 156 (1895) 

[29] H.C. Brinkman, "A calculation of the viscosity and the sedimentation constant for so- 
lutions of large chain molecules taking into account the hampered flow of the solvent 
through these molecules," Physica 13, 447(1947) 

[30] C.-O. Hwang, J. A. Given and M. Mascagni, "A New Fluid Permeability Estimation," 
to appear in Proceedings of the First Southern Symposium on Computing, (1999) 

[31] I.C. Kim, and S. Torquato, "Effective conductivity of suspensions of overlapping 
spheres," J. Appl. Phys. 71(6), 2727 (1992) 

[32] M. J. MacDonald, C.-F. Chu, P. P. Guilloit, and K. M. Ng, "A generalized Blake-Kozeny 
equation for multisized spherical particles," AIChE J. 37(10), 1583(1991) 

[33] L.M. Schwartz, N. Martys, D.P. Bentz, E.J. Garboczi and S. Torquato, "Cross-property 
Relations and Permeability Estimation in Model Porous Media," Phys. Rev. E 48, 4584 
(1993) 

[34] M. Doi, "A new variational approach to the diffusion and the flow problem in porous 

media," J. Phys. Soc. Jpn. 40, 567 (1976) 
[35] A. Cancelliere, C. Chang, E. Foti, D.H. Rothman and S. Succi, "The permeability of a 

random medium : comparison of simulation with theory," Phys. Fluids A 2(12), 2085 

(1990) 

[36] N.S. Martys, S. Torquato and D.P. Bentz, "Universal scaling of fluid permeability for 

sphere packings," Phys. Rev. E 50, 403 (1994) 
[37] R.E. Larson and J. J.L. Higdon, "A periodic grain consolidation model of porous media," 

Physics of Fluids Al(l), 38 (1989) 



17 



FIGURES 



2000 
1975 
1950 
1925 
1900 
1875 
1850 
1825 
1800 
1775 
1750 








0.05 



0.1 



0.15 



0.2 



0.25 



FIG. 1. CPU time required to calculate the capacitance of the unit cube, to a fixed tolerance, 
using the first-passage method. The variables 5h and 5i prescribe the distance from the surface of 
the cube that a diffusing particle must be, resp., to declare that a "hit", i.e., an absorption event 
has occurred, or to use a first-passage surface that intersects the surface of the cube. The figure 
shows that, for each value of Sh, there is a non-zero optimal value of 5i. The walk-on-spheres 
( WOS) method is given by 5i = 0. 
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FIG. 2. Permeability vs R/a for unit capacitance estimation methods applied to the randomly 
overlapping sphere model: the permeability estimate does not depend on sample radius even though 
with small sampling sizes more porous samples are needed. For R = 5, 10, 50 porous samples were 
used. When 4>i decreases, estimation is not stable due to the estimation instability. The solid lines 
are the average of the seven permeabilities from different sampling sizes. 
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FIG. 3. Permeability vs R/a for penetration depth estimation method applied to the randomly 
overlapping sphere model: the permeability estimate does not depend on sample radius, provided 
that this radius is chosen greater than a (porosity-dependent) minimum value. The minimum radius 
increases with porosity. 
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FIG. 4. Dimensionless permeability k/a 2 vs. porosity for the randomly nonoverlapping sphere 
bed (a = 1.0). The filled upper triangles are our simulation data with the effective radius method, 
the filled diamonds are our UC simulation data with the sharp boundary method, the filled circles 
are our PD simulation data with the sharp boundary method, the solid line is the Stokes' law, the 
long dashed line is the Happel-Brenner approximation, the squares are the Torquato-Kim upper 
bound, and the dot-dashed line is the Kozeny-Carmen relation. 
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FIG. 5. The dimensionless permeability k/a 2 vs porosity for randomly overlapping spheres of 
radius a = 1.0. The filled circles are our simulation data with the effective radius sampling method, 
the filled diamonds are our simulation data with the sharp boundary sampling method, the squares 
are the Torquato-Kim upper bound, the upper triangles are the Doi upper bound, the left triangles 
are the overlapping sphere bed data points from finite difference solution of the Stokes equation, 
and the long dashed line is the lattice-Boltzmann result. 
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FIG. 6. The dimensionless permeability k/a 2 vs porosity for SC spheres of radius a = 1.0 with 
R = 15.0. The diamonds and stars are our simulation data with the sharp boundary sampling 
method and the circles are L arson- Higdon calculation. 
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FIG. 7. The dimensionless permeability k/a 2 vs porosity for BCC spheres of radius a = 1.0 
with i? = 15.0. The diamonds and stars are our simulation data with the sharp boundary sampling 
method and the circles are Larson-Higdon calculation. 
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FIG. 8. The dimensionless permeability k/a 2 vs porosity for FCC spheres of radius a = 1.0 
with R = 15.0. The diamonds and stars are our simulation data with the sharp boundary sampling 
method and the circles are Larson-Higdon calculation. 
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FIG. 9. Dimensionless permeability k/a 2 vs porosity for SC lattice in the case of R = 50.0, 
a = 1.0. The stars are our penetration depth simulation data with the sharp boundary sampling 
method and the circles are Larson-Higdon calculations. 
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FIG. 10. Dimensionless permeability k/a 2 vs porosity for BCC lattice in the case of R = 50.0, 
a = 1.0. The stars are our penetration depth simulation data with the sharp boundary sampling 
method and the circles are Larson-Higdon calculations. 
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FIG. 11. Dimensionless permeability k/a 2 vs porosity for FCC lattice in the case of R = 50.0, 
a = 1.0. The stars are our penetration depth simulation data with the sharp boundary sampling 
method and the circles are Larson-Higdon calculations. 
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FIG. 12. Permeability vs R/a for penetration depth estimation method with the sharp boundary 
sampling method for SC lattice model. The permeability estimate does not depend on sample radius 
provided that this radius is chosen greater than a porosity-dependent minimum radius. 
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FIG. 13. Permeability k vs porosity for a matrix constituted by a polydisperse mixture of 
randomly overlapping impermeable spheres with the sharp boundary sampling method; the sphere 
radii are chosen to have the four values a = {1.5,3.5,5.5,7.5} with equal probability. Here the 
sample radius R = 50.0. 
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